Assessment of Higher-order RANS Closures in a 
Decelerated Planar Wall-bounded Turbulent Flow 

Elbert Jeyapaul* 

National Institute of Aerospace, Hampton, VA 23666 

Gary N. Coleman^ and Christopher L. Rumsey^ 

NASA Langley Research Center, Hampton, VA 23681 


A reference DNS database is presented, which includes third- and fourth-order moment 
bndgets for unstrained and strained planar channel flow. Existing RANS closure models 
for third- and fonrth-order terms are surveyed, and new model ideas are introduced. The 
various models are then compared with the DNS data term by term using a priori test- 
ing of the higher-order budgets of turbulence transport, velocity— pressure-gradient, and 
dissipation for both the nnstrained and strained databases. Generally, the models for the 
velocity— pressure-gradient terms are most in need of improvement. 
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Superscript 


+ 

Subscript 

(•) 


wall units, scaled by Ur and v 
normalization by the respective r.m.s. 

indices of tensors 
Reynolds-averaged quantity 


I. Introduction 

Second and higher-order Reynolds-averaged Navier-Stokes (RANS) closures are primarily motivated by 
the need to accurately predict complex flows. Convection-dominated flows, such as atmospheric boundary 
layers, have been shown to benefit from higher-order moment transport models;^ the third-order RANS model 
described in Hanjalic and Launder^ predicts the scalar flux accurately in a stably stratified flow, while the 
transport models used in second-moment closures (SMC) can capture the observed negative (‘up-gradient’) 
turbulent heat flux.^ The accurate prediction of the third-order moment transport model of Kurbatskii 
and Poroseva^ provides another example of the advantages of higher-order schemes. On the other hand, 
while higher-order closures are potentially more universal than their lower-order counterparts, they require 
more terms to be modeled and a larger set of coupled equations to be solved. The number of equations 
required increases rapidly as the dimensionality of the flow increases, as illustrated in figure 1. This cost is 
compounded by numerical stiffness"^ due to the difference in time-scales associated with moments of different 
orders. 



Shear Boundary Layer Boundary Layer 

Figure 1. Minimum number of transport equations necessary for higher-order incompressible RANS closures 
of different flows. In practice, additional transport equations are often used to account for the dissipation 
terms. 

In this paper we limit our attention to third- and fourth-order moments and the terms in their budgets 
(i.e., their transport equations), in an attempt to evaluate and improve the fidelity and universality of higher- 
order RANS closures. Beyond SMC, most of the research has been on truncation at third order, with an 
appeal to the Millionshtchikov hypothesis,^ which gives an algebraic closure for the fourth-order correlation 
as a product of quadratic moments. This hypothesis is exact for a Gaussian velocity distribution. As the bulk 
of most flows of interest are non-Gaussian, this approximation can be inaccurate. In the fourth-order-moment 
transport equations, models for the fourth-order velocity-pressure-gradient correlation and the dissipation 
are required, as are models for the fifth-order velocity correlations. In a recent hypothesis,® fifth-order 
moments have been approximated by the sum of products of quadratic and cubic moments. This hypothesis 
assumes the probability distribution function to be a Gram-Charlier series, i.e., a truncated approximation 
about the Gaussian distribution. The Gram-Charlier approximation provides good a priori predictions in 
boundary layers.® 
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For simplicity, we analyze wall-bounded turbulence in a parallel-flow geometry. This idealization allows 
us to consider fewer equations than the full spatial case (figure 1). For parallel flows, like this one, the 
third-order RANS closure requires the solution of five velocity-momentum transport equations while seven 
are involved in a fourth-order scheme. At least one additional scale-setting equation, such as for dissipation, 
is also used in practice. Reference data consisting of results from the direct numerical simulation (DNS) of 
channel flow before and after it is subjected to the straining field introduced by an adverse pressure gradient 
(APG)^’® are used to perform a priori testing of existing third- and fourth-order models. This work is 
thus in the spirit of that presented by Yorke and Coleman® and Sciberras and Coleman,^® who respectively 
performed a posteriori testing of scalar-eddy-viscosity and SMC schemes using this flow. Several new third- 
and fourth-order models are also introduced in this work. 

This paper is organized as follows. First, the strained-channel strategy used in the DNS is described, 
then the moment transport budget equations are provided. The DNS results are then summarized. Section 
IV includes both an overview and a discussion of higher-order modeling implications, where ‘higher-order’ 
refers to third- and fourth-order. A review of existing higher-order closure assumptions is then given, and 
several new modeling ideas are also presented. Finally, a priori model testing is conducted using the DNS 
results. By looking at the behavior of the various models in an a priori sense, we hope to gain some insight 
into which models might perform best when included in a complete model. It is recognized that a priori 
testing of model components has only limited practical value, because in many cases large terms almost 
balance each other and it is the combined effect of all the terms that matters. A posteriori testing of an 
entire scheme will ultimately be necessary to properly evaluate any chosen models. 

II. Strained-channel strategy 

Higher-order RANS closures are tested by comparing them with DNS of both conventional turbulent 
plane-channel flow and a strained-channel idealization of the APC boundary layer. The latter emulates the 
spatially developing low-Mach-number APC boundary layer by simultaneously applying in-plane ‘sliding’ 
motion of the walls and straining the domain of an incompressible turbulent channel flow, as shown in 
figure 2. The in-plane wall motion causes the bulk-flow deceleration needed to reduce the wall shear stress. 
The wall velocity Uw(t) is coupled with the mean streamwise centerline velocity Uc, such that their difference 
decreases in time, according to Uc — Uyj = Uc{0) exp(An), where Uc{0) is the mean centerline velocity at t = 0, 
when the strain is applied, and An < 0 is the streamwise compression induced by the (virtual) APC. (Recall 
that in a parallel flow, wall motion is physically equivalent to a spatially uniform acceleration/deceleration.) 
Consequently, in the frame of reference attached to the walls, the history of the mean centerline velocity is 
given by Cc(0) exp(Ant). 

The (rate of) applied strain Aij is irrotational and uniform in space; it steps from zero at t = 0 and is 
constant thereafter, which causes the entire domain (including the walls) to deform (figure 2). Here only the 
‘APC components’ are involved, namely the streamwise deceleration An < 0 and wall-normal divergence 
A 22 = —An > 0. For further details the reader is referred to Coleman et al.^ 



Figure 2. Turbulence in (a) a spatially developing APG boundary layer is simulated as (b) a time-evolving 
plane channel flow by combining the effects of irrotational strain that deforms the domain and sliding of the 
walls. 

The strained-channel approach has the advantage of producing the desired perturbation in an uncompli- 
cated parallel-flow geometry that contains essential features of the flow in question. Features of the APC 
boundary layer that are captured include the straining effect of the divergence of outer-layer streamlines and 
the reduction of wall shear stress. Missing is the effect of streamline curvature, which becomes important 
as the flow approaches the separation point. Also, pressure in an APC boundary layer experiences no ‘re- 
flection’ effects from the inviscid outer region, whereas in the channel the pressure and its correlations are 
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Figure 3. Spatially developing APG boundary layer represented by a temporally developing strained channel. 


influenced by the opposite wall. These limitations are offset by the fact that the present approach reduces 
the RANS problem from a two-dimensional (2D) steady flow to a one-dimensional unsteady analog, which 
allows efficient model testing. This strategy has been previously used to assess eddy-viscosity® and SMC^® 
models applied to the APG case. The reduction in dimensionality also greatly simplifies the higher-order 
RANS testing done here, since the parallel-flow assumption requires the solution of only half the number of 
velocity-moment equations required for the spatial counterpart (see figure 1). 

III. Moment transport budgets 

For the unsteady parallel-flow idealization of the APG boundary layer used here, for which the only 
nonzero mean spatial gradients are in the wall- normal direction (X 2 ), the second- moment transport equation 
reduces to 

dtUiUj = Vfj -\- + Tij + Tlij + Eij -\- T>ij, ( 1 ) 

where 

Vfj = —UiU2d2Uj — UjU^d2Ui 
Vfj = -WukAjk - uJu^Aik 

Tij = -d2U2UiUj 

Aij = -{1/ p){ujd^p + Uidjp) 

Eij = -2vdkUidkUj 
T>ij = i/d^UiUj, 

where 82 is the spatial derivative in the wall-normal direction. 
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( 2 ) 


The third-moment transport equation is given by 

dtUiUjUi = 'Pfji + 'Pfji + 'P^ji + Tiji + -|- Eiji + 

where 

Vfji = —UiUjU2d2Ul — UjUlU2d2Ui — UlUiU2d2Uj 
T^ijl — UjUlUf^A-il^ UlUiUk-^jk 

vfji = U^d2U^l + UjWld2U^i + UilHd2U2Uj 
= -d2U2UiUjUi 

liiji = -{l/p){uiUjdip + UjUidip + UlUidjp) 

Sijl — ‘2h'(^UidkUj0kUl UjdkUidkUl UldkUidkUj^ 
V^ji = vdluiUjUl, 


while the fourth-order moment counterpart is 


where 


dtUiUjUlUrn — 'P ijlm d” ^ ijlm d” ijlm d“ ^ijlm d- ^ijlm d- £ijlm d- 


'Pijlni = -UiUjUiU 2 d 2 Um ~ UjUiUmU2d2Ui - UiUjnUiU2d2U j - UmUiUjU2d2Ui 
^ ijlm ~ y>iUjUlUk-^mk UjUl^mlJ^k-^ik UlUmUiUk^j k UmUiUjUk-i^lk 

= UiUjUid2U2Um + UjUiUmd2W>H + Wh^AWU^j + UmUiUjd2W^l 

'Tijlm — ^2’1^2'l^i’l^j‘l^lUm 

Hijim = -{^/ p){UiUjUidmP + UjUlUmdiP + UlUmUidjP + UmUiUjdlp) 

£ijlm = -2v{uiUjdkUldkUm “f UiUidkUjdkUm + UjUldkUidkUm 

+ UjUmdkUldkUi + UlUmdkUjdkUi + UjnUidkUldkUj) 

'^ijlm — l^d2UiUjy’lUm' 


( 3 ) 


For future reference, we note that the velocity-pressure-gradient correlations n^-, 11^; and Ilijim can each 
be decomposed into a term {</>) involving the product of pressure and gradients of velocity and another (ip) 
containing only spatial gradients of pressure- velocity correlations, such that 11^ = (pij+ipij, 11^/ = (piji+ipiji 
and Uijim = 4 >ijim + i’ijim, where 


(pij = (1/ p)[pdiUj +pdjUi), (i.e. the pressure-strain term), (4) 

(piji = { 1 / p){pdiUiUj +pdiUjUi +pdjUiUi), ( 5 ) 

(pijlm = iM p){pdmUiUjUl + pdiUjUlUm + pdjUlUmUi + pOlUmU^Uj), ( 6 ) 

ipij = {—l/p){d2puj52i + d2pui52j), (i.e. the pressure transport term), (7) 

i’ljl = {-^/p){d2PUiUjS2l + d2PUjUlS2i + d2pUlUi62j), (8) 

Ipijlm = p){d 2 PUiUjUl 52 m + d2PUjUlUm^2i + d2pUlUmUiS2j + d2PUmUtU j 821) ■ ( 9 ) 


The terms on the right-hand side of (l)-(3) balance the total material derivative Dt. However, for the 
present parallel flows, the convection term is zero. With a constant streamwise pressure gradient and no 
applied strain, the statistics are steady, and Dt = 0. When the flow is strained and the pressure gradient is 
replaced by the effect of the in-plane wall motion , the material derivative is Dt = dt ^ 0, and the moments 
evolve in time, analogous to the spatial development of an adverse pressure gradient boundary layer. ^ In 
other words, the streamwise spatial derivatives of boundary layers are replaced here by time derivatives. 

In equations (l)-(3), terms on the right-hand side are referred to as (nominal) production due to shear 
7^'^, production due to Reynolds stress , turbulent transport T, velocity-pressure-gradient correlation 
n, dissipation e, and viscous diffusion D. The applied strain produces non-zero terms in all the moment 
equations, termed production due to applied strain . The terminologies used to describe the terms in the 
third- and fourth-order budgets are inherited from second-moment closures, and do not necessarily describe 
their physical behavior. That is, ‘production’ can be negative and ‘dissipation’ can be positive in the third- 
and fourth-order moment equations. 
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IV. DNS results 


A. Overview 

The DNS is performed by a pseudo-spectral (Fourier/Chebyshev-r) method. The parameters and procedures 
of this method are described in Coleman et al.^ The simulation was repeated for this study because the third- 
, fourth- and fifth-order moments of interest here were not part of the original study. The new statistics were 
evaluated from a new ensemble of 159 fields (in each field, averaging was also performed over streamwise- 
spanwise planes and folded about the centerline, taking advantage of the symmetry), using 256 x 193 x 192 
spectral modes in a 27 ti5(0) x S{0) x tt6{0) domain (where (5(0) is the half-height of the initial/pre-strained 
channel), in the streamwise x = Xi, wall-normal y = X 2 and spanwise z = x^ directions, respectively. The 
runs were made on Pleaides, an SGI® ICE-X cluster using Intel Xeon processors at the NASA Advanced 
Supercomputing division. The unstrained channel simulations required 11,000 processing core hours to 
obtain the 159 independent realizations, which were used as initial conditions for the corresponding 159 
strained-channel runs, requiring a further 8,000 (« 50 x 159) hours to advance each of them to a strain time 
of A 22 t = 0.77 (well past the point of separation; see below). 

The fully-developed (unstrained) channel Reynolds number i?e,- = Ut5/v is 392, which is large enough 
to sustain a well-defined inertial sublayer. At this initial unstrained state (^ 22 ^ = 0), the Reynolds number 
based on the mean centerline velocity Rcc = Uc^jv is 7,910, while the momentum-thickness and bulk 
Reynolds numbers are respectively Re% = Uc9lv=7Qi and Rcm = 25C/m/i^=13,770 (where 9 is the half- 
channel momentum thickness, and Um the bulk, wall-to-wall, average velocity). As in Ref. 7, a relatively 
weak APG is applied to the series of realizations of unstrained channel-flow turbulence, with A 22 = —^11 
chosen to be 31% of Mt( 0)/5(0), the ratio of the initial friction velocity to the initial channel half- width. The 
applied strain |An| is at least an order of magnitude smaller than the mean shear dU/dy over the entire 
channel for all strain times.® Consequently, the magnitude of the applied irrotational strain is not large 
enough at any channel height to cause the non-linear interactions of turbulence to be negligible; it is these 
nonlinear interactions that the higher-order RANS models are designed to predict with higher fidelity than 
SMCs. Both the strain-dependent ‘fast’ and strain-independent ‘slow’ parts of the budget terms are thus 
expected to contribute to the moment evolution. 

The resulting mean velocity profiles are shown in figure 3. Many of the basic features of APG boundary 
layers are present. For example, the shape factor (in terms of the displacement and momentum thickness 
in the half-channel) increases from H = S* /9=lAb at ^ 22^=0 to H = 1.7 at ^22^=0. 365 and H = 2.5 at 
^ 22 ^= 0 . 77 - past the time at which the mean ‘separation’ (i.e. mean-flow reversal at the wall) occurs, at 
A 22 t = 0.68. This is consistent with an axisymmetric-body separation-bubble,’ for which H « 2.7. The 
effective Clauser parameter j3eff = —5*UcAxi/ul is 0.78 at ^ 22 ^ = 0, before increasing to 5.7 at ^22^=0. 365. 

B. Higher-order-modeling implications 

Normalized statistics such as skewness and flatness of velocity correlations vary through the channel height 
(figure 4). Higher-order models will need to resolve this variation by using models for unknown terms to 
close equations (2) and (3). The unstrained-channel skewness and flatness profiles compare well with prior 
channel-flow DNS.^® Experimental studies of a turbulent boundary layer^^ show the skewness and flatness 
reach their global maxima or minima near the wall and the boundary-layer edge, with a nearly constant 
value between the two extrema. This trend is not noticed in the channel- flow skewness (figure 4a), in that 
here both S\ and S 2 do not have their maxima/minima near the wall nor at the channel centerline. In the 
inertial subregion, the APG strain (until ^ 22 ^ = 0.365) drives the skewness for all three velocity components 
toward zero (i.e. toward the Gaussian state). The strain also causes the flatness for the wall-normal (F 2 ) 
and spanwise (F 3 ) velocities to reduce throughout the channel toward the Gaussian value of 3, and for the 
streamwise component to become slightly less Gaussian than it was before the strain is applied, reducing 
from Fi « 3 to about 2.7. The effect of the straining in this region is similar to that seen for the mild- APG 
boundary layer^^ (the skewness and flatness differences observed in the very-near- wall and wake regions are 
presumably related to near- wall measurement difficulties and/or to the difference between the channel and 
boundary-layer geometries) . 

For the present flow, the moments that act as dependent variables in a fourth-order RANS closure 
are —uv, mu®, m^, uv^ and (Models for the pressure- velocity correlations and dissipation could 
in principal introduce dependence on additional moments, but we do not consider that possibility here.) 
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Figure 4. Plots of (a) skewness Sa = and (b) flatness Fa = u'^Ku'^)'^ (no sum over a = 1,2,3): , 

A22t=0; , j422l=0.365; o, ui = u; V , U 2 = v and A, ug = w. 


The evolution of the relevant third- and fourth-order moments and the terms in their budget are shown in 
figures 5-8 for the unstrained case ^22^ = 0 and for the strained case at 2)22^ = 0 . 365 . The smooth nature of 
the profiles in figure 4 indicates that the statistical error is fairly low. For the unstrained flow, checks were 
performed to ensure that these statistical errors are negligible, by comfirming that the sum of the right-hand- 
side budget terms from the 159 -field ensemble is much smaller than the smallest term in each budget. It was 
also verified that the third- and fourth-order moments were not corrupted by aliasing errors, by recomputing 
them after projecting the velocity fields onto a collocation grid that is 2 ^ times larger, and confirming that 
no significant differences result, compared to the moments obtained from forming the products on the grid 
used for the DNS. Moreover, the sum of the right-hand-side terms in equations ( 2 ) and ( 3 ) are very close to 
the observed time rate of change of the moment being transported (figure not shown), which gives further 
confidence in the DNS and the calculations of the terms in each of the higher-moment budgets. 

The third-order moments, which measure the skewness of the corresponding probability distribution 
function (PDF), are smaller than the second- and fourth-order moments. The third-order moments and 
—uv'^ have inflectional profiles, with negative values near the wall and positive values elsewhere (figures 5 and 
6 respectively). The left-hand-side subplot in each of these composite figures shows the moment evolution 
and, on the right-hand side, the strain- induced changes to each of the terms in its budget, equation ( 2 ). We 
notice negative values at yw /3 < 0 . 1 , which is primarily caused by an Imbalance between II222, 'P222 
7222 - Both n and T are unclosed terms, needing modeling. Note that V§22 = 0 ; hence is not affected 
directly by the mean velocity. The direct effects of the applied irrotational straining are manifested through 
T’^2- The strain applied at ^122^ = 0 causes an instantaneous increase in both the velocity-pressure-gradient 
term (due to incompressibility) and the additional negative ‘production’ As seen in inset figure 5 d, 

■p^2 dominates II222J causing to increase throughout the channel instantaneously. As noted in Coleman 
et al.,' this sudden increase is due to the channel turbulence not receiving a ‘warning’ of the impending 
discontinuous temporal change. This behavior is thus unlike a spatially developing boundary layer where 
convective changes In the mean flow propagate upstream. However, at ^22^ = 0 . 365 , the difference between 
the inner- and outer-layer dynamics causes to decrease in the inner layer and grow in the outer layer, 
similar to the Reynolds stresses.^ Notice (compare inset (c) to inset (d) of figure 5 ) that dv^/dt is of the 
same order of magnitude at ^22^ = 0 and 0 . 365 . This implies that the step increase of An = A22 < 0 
at 0 does not derail the turbulence mechanisms, allowing the nonlinearities (manifested through change in 
budget terms) to take effect at later times, causing the different inner-outer layer dynamics. The inner-layer 
reduction and outer-layer growth of are also seen past flow reversal (^22^ = 0 . 77 ). 

The evolution of uv"^ is caused by active contributions from nearly all the budget terms (figure 6 ). An 
inner-layer decrease and outer-layer growth is noticed, similar to that for v^. The direct effect of straining 
through V122 has a lesser role in duv^ jdt than the net contribution (indirectly affected by straining) from 
the rest of the budget terms. 

The fourth-order moments, which are all positive definite, show noticeably larger outer-layer increases 
than inner-layer decrease (figures 7 and 8 ). This difference between inner- and outer-layer response is thus 
different from that for the third-order moments. Interestingly, we notice responds to strain in a manner 
similar to ^ although an additional term ( 7 ^^) contributes. The v‘^ budget is driven primarily by an 
imbalance between H2222, £2222 and T2222 (figure 7 ). The same budget terms contribute to the evolution of 

(figure 9 in Ref. 7 ). The wall- normal flatness F2 increases rapidly near the wall (figure 4 ) due to the 
differences in the rates of decay of the dominant budget terms, H, e and T, of which H2222 and £2222 remain 
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Figure 5. (a) Third— order wall— normal velocity correlation (r®) profiles: , A 22 t = 0; , A 22 t = 


0.365; , A 22 t = 0.77. (b) Terms in the budget at A 22 t = 0.365: , Production due to shear 


(pS). ^ 

Production due to Reynolds stress 



5 

Turbulent transport 

(T); , 

Velocity pressure— gradient 

(n); , 

Dissipation 

(e); -■■■ 

, Production due to applied strain (P^), 


also shown in insets (c) and (d); , Viscous diffusion (T>). Thin curves denote terms at t=0 (before 

strain) and are identified by the shaded regions, which indicate change from unstrained initial conditions. 
The bndget terms in (b) are normalized by u^(0)/u. Inset (c) Hollow circles (O) are snm of all budget terms 
(Ri dv^/dt) at A 22 t = 0.365. (d) Budget term balance dne to instantaneous application of strain {A 22 t = O'*"). 
Curves shown are snbtracted from initial nnstrained profile (A22t=0). The curve and symbol types are as in 
(b) and (c). The axes in insets (d) and (c) are non-dimensionalized as (b). 






Figure 6. (a) Third— order shear— stress correlation (— nr^) profiles: , ^ 22 ! = 0; , A 22 t = 0.365; , 

A 22 t = 0.77. (b) Terms in the budget at A 22 ^ = 0.365. Inset (c) Hollow circles (O) are sum of all budget 
terms (ps duv^/dt) at A 22 ^ = 0.365. Inset (d) budget term balance due to instantaneous application of strain 
(A 22 t = 0"^). The budget terms are normalized by u^{0)/i'. For line legend and shading key refer to figure 5. 
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Figure 7. (a) Fourth— order wall-normal velocity correlation (v^) profiles: , A 22 t = 0; , A 22 t = 0.365; 

, A 22 t ~ 0.77. (b) Terms in the budget at A 22 ^ ~ 0.365. Inset (c) Hollow circles (O) £^re sum of all 

budget terms (ps dv^/dt) at A 22 ^ = 0.365. Inset (d) budget term balance due to instantaneous application of 
strain {A 22 t = O'*"). The budget terms are normalized by u^{0)/u. For line legend and shading key refer to 
figure 5. 



Figure 8. (a) Fourth— order shear— stress correlation {uv^) profiles: , ^ 22 ^ = 0; , A 22 ^ = 0.365; , 

^ 22 ^ = 0.77. (b) Terms in the uv^ budget at A 22 f = 0.365. Inset (c) Hollow circles (O) are sum of all budget 
terms duv^/dt) at ^ 22 ^ = 0.365. Inset (d) budget term balance due to instantaneous application of strain 
(A 22 t = 0"^). The budget terms are normalized by u^{0)/iy. For line legend and shading key refer to figure 5. 
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the largest, since they decay slowest as > 0, being 0{y%) (see Appendix). It is therefore important 
that models for II 2222 and £2222 reproduce this near-wall asymptotic behavior. The direct effect of straining 
through the applied production (7^^22 = —4:VvvvA22) nearly balances the net effect of the other three 
dominant terms mentioned earlier (figure 7), causing v'^ to grow for y^/S > 0.3 and decay near the wall. 
The centerline v‘^ decreases between ^ 22 ^ = 0 and 0.365 and increases from ^ 22 ^ = 0.365 to 0.77, as the flow 
slowly recovers from the abrupt initial strain. 

From the uv^ budgets in figure 8 we notice all the budget terms make a significant contribution to 
duv^/dt. Production due to mean shear is a dominant term in the balance. Production due to turbulence 
(^ 1222 ); although small, couples the uv^ evolution to the third-order moments (a route for non-Gaussian 
effects to manifest themselves). In the near-wall region {y^/S < 0.1) the evolution of uv^ is primarily driven 
by ^1222 • 

Since Rotta,^^ it has become conventional in SMC to split the velocity-pressure-gradient term (11^^) into 
the sum of pressure-strain 4>ij and pressure-transport (often referred to as ‘pressure diffusion’) ijjij terms (see 
earlier equations (4) and (7)). For modelers this is advantageous, since pressure-strain is redistributive/trace- 
free (i.e., —(I)ii = (j )22 + 4 > 33 ) j and is thus responsible for transfer of energy between the three components, and 
pressure-transport is non-zero only near the wall. Therefore, in homogeneous turbulence, Uij = (f>ij, making 
the problem amenable to the use of homogeneous turbulence theory to develop models. By definition, in 
homogeneous flows, ipiji = ipijim = 0, hence homogeneous theory can be extended to model the higher-order 
counterparts of On the other hand, the wall boundary conditions on (jjij and are both non-zero, 
in fact they are equal and opposite, and therefore more complicated than the zero condition on each of 
the elements of n^-. This difficulty is a symptom of the larger problem that in the near- wall region, the 
magnitudes of (f>ij and ipij are much larger than their net effect, which will tend to amplify any errors in 
their respective models. 

Some modelers^ A4 j^ave applied the II splitting to third-order budgets. In contrast to the second-order 
case, (p = 0 and = 0 at = 0 for both third and fourth order, and the fourth-order ‘pressure-strain’ is 
not redistributive, in the sense that pijij = {2/ p){pdjUiUjUj + pdiUiUjUj) ^ 0. Furthermore, as seen in the 
DNS results in figures 9 and 10, the magnitudes of (pm and ipm, and of (pijim and 'rpijim, are comparable 
in regions well away from the wall (although the APG strain does tend to drive tpijim ~ 0 there); we note 
that of the components considered here, only '02222 is nearly zero and II 2222 ~ 02222 in the outer-layer. And 
perhaps most importantly, from a modeling point of view, the wall-normal variations of 11^;/ and Ilij/m are 
simpler then that of their constituent 0 and 0. All of this suggests there is little if any inherent advantage 
to the n = 0 -|- 0 decomposition for the higher-order moments. 






fc/<5(i) 


Figure 9. Third— order velocity— pressure-gradient budget term split before and after straining. The velocity— 
pressure-gradient term II is ; pressure-strain term (f) is ; pressure transport term ^ is . 

Models for fourth-order moment evolution have been the least explored. Single-point fourth-order mo- 
ments have major symmetries {uiUjUiUm = UjUiUmUi = uiUmUiUj = UmUiUjUi) by definition. This reduces 
the number of independent components from 81 (= 3^) to 21 and the number of principal invariants from 
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y»/<5(t) 



y»/<5(i) 


Figure 10. Fourth— order velocity— pressure-gradient budget term split before and after straining. For the line 
legend refer to figure 9. 


9 to 6 (see Batten^®). The invariants are objective properties of the tensor, as they do not change with 
coordinate frame rotation. The simplest of invariants, the first, is the trace of the tensor, which for the 
fourth-order moment is tr(uiUjUiUm) = UiUjUiUj = u‘^ + v'^ + w'^ + 2{u'^v^ + v'^w^ + w^u^) . Turbulent kinetic 
energy is half the trace of the Reynolds stress tensor K 2 = {l/2)uiUi. Similarly, a fourth-order analog of 
turbulent kinetic energy can be defined as K 4 = {l/ 2 )uiUjUiUj. 

The quantity K 4 is positive definite and is a measure of the magnitude of the fourth-order moment. 
Its budget terms are descriptive of the terms they represent, that is, ‘production’ is positive definite and 
‘dissipation’ is negative definite (figure 11). We notice from figure 11(a) that as a result of straining, K 4 is 
reduced drastically in the inner layer and increases in the outer layer of the channel, similar to turbulent 
kinetic energy K 2 (figure 7 of Ref. 7). As for K 2 , for which the balance of budget terms is Sk ~ as 
yw — >■ 0, the dissipation versus viscous-diffusion balance is also observed for K 4 . However, in the K 2 budget, 
the dissipation and viscous diffusion are non-zero at the wall. We also notice that dtK 4 evolves similar to 
dtK 2 (compare figure 11(d) to figure 7 of Ref. 7). Evolution of K 4 in the outer half-channel is driven by 
applied-strain production (= 4(u"‘ -\- v'^ + vP-w"^ — v‘^w‘^)A 22 ), causing K 4 to steadily increase in the outer 
layer (figure 11(a)). The equivalent of a second-order turbulence time scale (t 2 = 2 K 2 /\eu\), a measure of 
the rate at which the turbulence kinetic energy is dissipated at the smallest scales, can also be defined for 
the fourth-order moment as T4 = 2AT4/| A plot of the ratios of these time scales (figure 11(c)) shows 

that T 2 ~ 2 t 4 for both the unstrained and strained cases throughout most of the channel, except near the 
wall. That is, fourth-order moments tend to respond on the order of twice as fast as second-order moments. 
Scale-similarity models (see below) for the fourth-order budget terms need to account for this more rapid 
evolution. A simple model relating the fourth-order timescale to second-order timescale can be given by 


T2 


T 4 = f^=2 + A exp 

Jw 



( 10 ) 


where, y'^{t) = (0)5(0) /5(f) is the wall-distance scaled by the initial unstrained channel height. The 

timescale ratio f-oj is only a function of wall distance, whose variation has been chosen to approximately 
reproduce the near- wall behavior of the DNS. 


V. Closure assumptions 

In this section, we introduce currently available models utilized in third- and fourth-moment transport 
RANS schemes. These modeling assumptions are summarized in Table 1. We also introduce several new 
models (using the nomenclature JCRl-6) in this section. Most of the models in Table 1, as well as the new 
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Figure 11. (a) Fourth-order turbulent kinetic energy analog {K 4 = UiUjUiUj) profiles at strain times: , 

A 22 t = 0; , A 22 t = 0.365; , A 22 t = 0.77. (b) Terms in the K 4 budget at A 22 t = 0.365. The budget 

terms are normalized by u^(0)/u. For line legend and shading key refer to figure 5. (c) Ratio of second- and 

fourth-order turbulent time scales (t 2 /'T 4 ), where T 4 = UiUjUiUj / at strain times: , A 22 t = 0 ; , 

A 22 ^ — 0.365. 


JCR models, will be compared to the DNS data in Section VI. 


Table 1. Existing models for terms in higher-order RANS. 


Moment-order 

r 

n 

5 

3 rd 

M5, KSKII 6 , CCH,1^ GHRSi®, 
GNi9 

DL1,3 NT, 14 KP1,4 KSK216 

Ai, DL2,33, KSK 3 I 6 

4th 

JDJ® 


KP220 


A. Turbulence transport 

One of the motivations of higher-order RANS is the ability to explicitly account for the turbulence transport 
term by solving an evolution equation of a higher-order moment (equations (2) and (3)). In contrast to 
SMC, for which the triple-velocity correlation terms have been approximated by various approaches - with 
limited success^ - for third- and fourth-order transport schemes, the task of modeling turbulence transport 
is delegated to approximating the fourth- and fifth-order moments. The models for T listed in Table 1 can 
be classified as PDF approximations, gradient-diffusion models, scale-similarity models (universal profiles 
in proportion to a relevant time scale), or a combination of these. The PDF and scale-similarity models are 
local models, as they relate the higher-order to lower-order moments. Unlike in second-moment equations, 
the turbulence transport term plays a dominant role in third-moment transport equations. In the near-wall 
region its contribution is of leading order (figures 5-6) . The transport terms involve gradients of fourth-order 
moments. Hence, high-fidelity fourth-order moment models are essential to predicting the evolution of triple 
correlations accurately. 

The widely used fourth-order moment model by Millionshtchikov,® which makes use of a quasi-normal 
(QN) hypothesis, assumes the velocity fluctuations have a Gaussian distribution. It relates the fourth-order 
moment to second-order moments: 

M: UiUjUlUm = UiUjUlUmlQN = WUj ■ UlUm + UjUl ’ UiUm + WM ‘ UmUj. (II) 

(The label at the start of the above, and following, equations corresponds to the entry in Table 1.) Equa- 
tion (11) is independent of skewness, as odd correlations are zero for a Gaussian variable. Implicit in this 
model is the assumption that the flatness is three. 

Gorrections to M have been proposed by a number of researchers, including modeling the cumulants 
(measures of the deviation from Gaussianity u'^ — m^Iqjv) as a transport process. An algebraic model for the 
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cumulant transport by Kawamura et al.^® approximates the cumulant as a gradient-diffusion process: 


KSKl: UiUjUiUk - UiUjUiUklQN = CrUkUmT-^ — u^UjUi, (12) 

OXm 

where Ct = 0.02 has been evaluated using the unstrained channel DNS. The turbulence time scale (r) is 
defined as t = k/e. Unlike the even-order moments, triple correlations have inflectional profiles (figures 5, 
6). Hence their gradients greatly influence the cumulants. 

The non-local model of Cheng et al.^^ for wall-normal fluctuation accounts for the distribution of triple 
correlations throughout the domain: 


CCH: 



= v"^\qn + Cl 



(13) 


where Ci = 4 and y = yw/^- Triple correlations are much smaller in magnitude than fourth-order correla- 
tions, and are not positive-definite. Hence corrections through this model are expected to be insignificant. 
This model is easy to implement because the turbulence transport term is the gradient of the fourth-order 
correlation. 

The local model by Gryanik et al.^® accounts for the non-Gaussian nature of the flow by making the 
fourth-order moment a function of the skewness. The model is developed as a solution to the 16-delta PDF 
function modeU® that satisfies additional constraints of tensor invariance, symmetry and realizability of the 
flow variables. The modeled moment is UiUjUiUm = a UiUjUiUmlQN, where the coefficient a is unique for 
each component, given as 


GHRS: 


— 1 + gS'?; a:^ — 1 -I- 


1-k 


Cuv 


-^1^2, 


(14) 


where Sa is the skewness of Ua, Cuv = uv, and the circumflex ‘ * ’ denotes normalization by the respective 
root-mean-square (RMS) value. The cumulant transport here is linear in UiUjUiUmlqN, leading to a large 
correction in regions where the product of QN flatness and skewness is high. Also, the correction is only 
positive for autocorrelations as it involves the square of skewness. 

Grossman and Narayan (1993)^® proposed a model for the wall-normal component, similar to GHRS, 
with 

GN: = 0.767 -f0.6S'|. (15) 

This model produces flatness F 2 < 3, contrary to what is observed from boundary-layer DNS® as well as our 
current channel DNS (figure 4). 

Realizability constraints for third-order correlations can be found from the Cauchy-Schwarz inequality, 
which states that {ujUjY ^ 


{uiUjUkY ^ min ' 


2 , 

(22 

2 

2 

ut 1 




2 I 

( 2 2 

2 

2 


[u^Uk 


■«fc 

< 

Uu] 


■“1 


)■ 


(16) 


This constraint relates the flatness {F) to the skewness (S'), giving F > 1 + S^. It can also be used to relate 
the third- and second-order moment by invoking the QN hypothesis. Realizability limiters can be designed 
for algebraic or differential-transport models of higher-order moments. 

The turbulent transport terms in the fourth-moment transport equations have a less significant role than 
they do in the third-order moment transport (cf. figures 5 and 6 with 7 and 8). The transport terms involve 
gradients of fifth-order moments. Models for the fifth-order and higher moments have been proposed by 
Jovanovic et ah,® in terms of a Gram-Charlier series-expansion of various orders for the probability density 
distributions of the turbulence fluctuations about a Gaussian/normal PDF (with the accuracy presumably 
increasing with increasing order of truncation). The fourth-order truncation is given by 
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u'^v = Qv?v ■ v? + 4uv ■ 


y3y2 — Qy^y . y2 y _|_ ^y2 . yy2 y2 . y3 ^ 

y2y3 — Qyy . yy2 _|_ 0^2 . y2 y _^y2.y3 

JDJ: (17) 

uv^ = + 4 mv ■ v^, 

= lOu^ • v^, 

= lOu^ • u^. 

The fifth-order moments from the DNS are inflectional.® Note that equation (17) assumes the inflectional 
points of the fifth-order moments to be same as those of the third-order moments. This assumption agrees 
reasonably well with DNS, as we notice second-order turbulence transport term (figure 9 in Ref. 7) to have a 
similar trend as fourth-order turbulence transport term (figure 7). The non-dimensionalized values of fifth- 
order moments are represented as super-skewness, SS = ft®. The magnitude of the fifth-order skewness is 
larger than the third-order skewness (|SS| > |S|) for all velocity components. All the terms in equation (17) 
have correct near-wall asymptotic behavior as — >■ 0. The JDJ model has been shown to provide a good 

approximation for the higher-order moments in a channel flow® and in a longitudinally rotating pipe.^® 
Nagano and Tagawa^^ have developed a ‘structural’ model for the triple correlations based on the sweep 
and ejection mechanisms (i.e. the so-called Q4 and Q2 events, respectively) in near-wall turbulence. The 
model relates the third-order ‘shear’ stresses to the skewness as 


yniyU = (J 


+ <yn¥ 


where {n,m) = (1,2) or (2,1) and 

C=- 


(f) -1 


with 


— a; is even, 
1, X is odd. 


(18) 


Hence, while modeling third-order moments, only transport equations for u® and u® need to be solved and 
triple cross-correlations such as uv"^ and v?v can be predicted from equation (18). 

Based on this idea, here we define a new model for turbulence transport terms in the fourth-order moment 
equations. Since the fifth-order correlations are structurally similar to the third-order correlations (as they 
measure the imbalance between Q4 and Q2 contributions), Nagano and Tagawa’s structural model (18) can 
be extended to write the fifth-order cross-component moments in terms of the super-skewness as 


JCRl: 


ymyn = (J 


amU^ + 


(19) 


where m-|-n = 5, and the model coefficients C, am and (j„ are the same as in equation (18). This new model 
can be used to predict the cross-correlations using autocorrelations, by either solving for transport equations 
for and u® or by using an algebraic model that relates super-skewness to skewness. We will include this 
new JCRl model in the comparisons in Section VI. 


B. Velocity pressure-gradient correlation 

As mentioned earlier, a few modelers extend the splitting of velocity-pressure-gradient (into pressure-strain 
plus pressure transport) to higher-order moments. Pressure-strain is modeled by Dekeyser and Launder® by 
the same strategy typically used in SMC. There is a linear return-to-isotropy term (‘slow’ part) and a set of 
production terms (‘fast’ part): 


DLl: 


= -Ci 


UiUjUl 




( 20 ) 


with constants C\ = 13.33 and C 2 = 0.5. Nagano and Tagawa^^ proposed another model for pressure-strain. 


NT: _|_ Cu2TUiUm 


a 


ul ^ 


dXr, 


1 


- -^k 


(21) 
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with the model constants Cui = 0.12 and Cu 2 = 1-4. The second part in this model is not directly dependent 
on strain (there is no ‘fast’ part). Hence response to changes in mean strain are indirectly resolved through 
changes in Reynolds stress, potentially lagged. Moreover, the second part is identical for the <(>iii, 4'ii2, 
and </>i 22 components, indicating these components have similar dynamics. The model of Kurbatskii and 
Poroseva^ for the velocity-pressure-gradient of the third-order moment is given by 

KPl: = b^ji 0>rnijln^nUjn- 

This model consists of a ‘slow’ part (biji) and a strain-dependent ‘fast’ term. The ‘slow’ part is a function 
of three coefficients and the tensor amijin is dependent on two model coefficients. This model will not be 
tested as the full expression of the model is not available in the open literature. 

An improvement to the model of Dekeyser and Launder^ was proposed by Kawamura et al.,^® in terms 
of the velocity-pressure-gradient: 

KSK2: (1 - /^2) + /^2, (23) 

where 

T 

with empirical constants Ci=0.3, Ci=0.6 and C2=0.1. In equation (23), a modified (piji was provided 
(compare with equation (20)), and a correct wall-asymptotic behavior was included by damping components 
as a function ( fw 2 ) of v^. The model assumes H^-/ = (piji away from the wall and H^-; is only dependent on 
a return-to-isotropy term in the near-wall region. 

An approach to modeling the velocity-pressure-gradient correlation is proposed here. While not a model 
per se (because it relies on an additional unspecified model for Tiji), this approach will be tested in an a 
priori sense in Section VI, with Tip taken from the DNS. The approach is derived from equation (2), and 
models the H^-/ term as 

JCR2: = Ca - a + Tp) . (24) 

The first term on the right-hand side represents the ‘slow’ term, where Cq/t accounts for the turbulence 
response, and the remaining terms account for the ‘fast’ part. The ‘fast’ part only has contributions from 
production due to shear and turbulence and turbulence transport, as they are dominant (figures 5 and 6). 
While dissipation has a role in the budget transport, its model expression (discussed in the next section) is 
not well known. Therefore, its effects are resolved instead through tuning the coefficients, chosen as Ca = 0.5 
and Cb = 0.9. 

Currently, there exists no model for the fourth-order velocity-pressure-gradient tensor. We propose an 
approach similar to the description in equation (24), with 

JCR3,4: n,p^ = Ca - Cb {Vfp^ + Vjp^ + %p„r) . (25) 

As with JCR2, this is not a complete model because Tijim is currently taken from the DNS. The effects of 
diffusion are again accounted for by the coefficient Ca- This choice is based on the observation that dissipation 
is self-similar to the moment it transports (refer to the next section). This approach was tested for fourth- 
order correlations, and we find that including the ‘fast’ part improves the predictions for wall-normal auto 
and cross-correlations. The simpler approach, JCR3, only accounts for the ‘slow’ part in equation (25) by 
using Ca = 2, Cb = 0, that is, assuming similarity with the fourth-order moments. The JCR4 version 
accounts for both the ‘slow’ and strain-dependent ‘fast’ part with Ca = 1.5, Cb = 0.3. Note that Ca < Cb 
in equation (24), while Ca > Cb in equation (25), indicating that IVipm is more scale-similar to the moment 
it transports than is Hip. 

C. Dissipation 

The majority of the models for dissipation in SMC assume isotropy and apply a near-wall correction to 
approach the wall values. One advantage to modeling of high-order dissipation is that the wall values of 
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dissipation for moments greater than second-order are zero. This leads to a simple boundary condition for 
dissipation. However, the near-wall asymptote as > 0 needs to be modeled through damping or other 
means. Assuming local isotropy of third-order moment dissipation leads to £^7=0; such a model was used by 
Andre et al.^ (A). From the DNS, we notice that such a model is grossly incorrect, so it will not be included 
in the testing below. 

A more accurate model, by Dekeyser and Launder,^ relates Siji to its lower-order term by invoking a 
generalized gradient-diffusion hypothesis^^ (GGDH). They further assume isotropy of second-order moment 
dissipation, that is Stj = {2/3)sSij, with e = \eu\ to relate it to the scalar diffusion. The model is 

ds 

OX]i 

with Cs = 1. However, equation (26) results in Sm = 3ei22 and £222 = 3£n2, whereas the channel DNS 
shows these values to be unequal in magnitude and sign. Note that third-order moment ‘dissipation’ is not 
dissipative, as the components are not monotonically negative in the channel flow.® The inflectional profile of 
Eiji will not be reproduced by this model because the gradient of isotropic dissipation is always non-negative 
in a channel flow. This shortcoming can be overcome by reverting to the GGDH model, without assumptions 
of isotropy on the second-order dissipation. Such an anisotropic model of DL2 (with Cg = 1) is 

T-.TO. f , 9e,k , dejk\ 

DLo. ^ijl — I “t“ UiV-l J ■ (2' ) 

A model based on similarity of Siji to the third-order moments was developed by Kawamura et al.^® The 
expected near-wall damping (refer to the Appendix) is provided using the invariants of the Reynolds stress 
anisotropy and the componential damping using a wall-normal vector (rii). The model (with constant Cg = 2) 
is 

KSK3: £„, = (1 - /^i), (28) 

where, 

e*jl = - {iUiUjUl + 2 {u^UjUmrimni + UjUlUmnm.ni + UlUiUmnmrij) 


The GGDH has been extended to the fourth-order dissipation tensor by Kurbatskii and Poroseva.^° The 
resulting expression for dissipation is 


KP2. Eijlm — C JE4: 


f dEi^, 


dEji 

\UiUjUk 

-\-UiUlUk ^ 

OXk 

+ UrUmUk^ 

OXk 


Oeu 

dEij 

+ Uj ai Ilk p. 

OXk 

^m'^k 0 

OXk 

H" ^k 

OXk 


(29) 


We find Cea to not be a constant for the components of dissipation; however, a value of 2 gave a reasonable 
fit to the DNS (refer to following section). 

In the fourth-order budgets, the dissipation term Eijim is always dissipative (figures 7 and 8). We 
propose a new model based on similarity of dissipation to the fourth-order moment. One of the weaknesses 
of such a model is that free-stream values of dissipation are incorrect, as at the channel centerline UiUjUiUm 
are non-zero while dissipation Eijim is zero. Two variants are proposed, where the fourth-order timescale 
is T 4 = UiUjUiUj / Eimim and T 2 is based on the second-order timescale (t 2 = utUi/Ejj). The constants 
Cl = 2 = 2 C 2 are recommended: 


JGR5 : 

^ijlm 

UiUjUlUm 
= Cl — ^ 


T2 

JCR6 : 

^ijlm 

UiUjUlUm 

= C 2 — 


T4 
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VI. A priori model testing 


In this section, closure- mo del components of all u and v correlations of third- and fourth-order moments 
are tested a priori - that is, the individual closure assumptions are compared to the corresponding DNS 
results by using component DNS data as input to the models. The limitations of this approach are ac- 
knowledged, namely that a RANS scheme’s overall behavior is due to the combined effect of all the closure 
assumptions, which in general will include many compensating errors, such that a scheme’s net prediction 
may be much better in practice than implied by the validity of the individual assumptions. A complete 
picture will require a posteriori tests of an entire scheme (as in Sciberras and Coleman^®). This limita- 
tion notwithstanding, a model developer will presumably benefit by knowing the accuracy of the individual 
closure assumptions, which motivates the present study. The main goals are to assess: (i) the ability of 
the individual closure assumptions to reproduce the statistics of the DNS, (ii) the effectiveness of model 
coefficients, tuned for an optimal fit of all components,^® and (iii) the merits of general modeling concepts, 
and thereby their potential usefulness for 2D attached and APG flows. 

A. Unstrained channel 

The unstrained channel flow is a steady fully developed channel flow. That is, the moments of the statistics 
do not evolve with time, as the budget terms balance each other, leading to 5* = 0 in equations (2) and (3). 
The flow is in turbulence ‘equilibrium’, with the production of turbulence kinetic energy roughly balancing 
dissipation in the (nominal) log layer of the flow.^® We also notice ‘equilibrium’ in the fourth-order turbulence 
kinetic energy analog (hgure 11). The statistics presented here are representative of a moderate Reynolds 
number {Rct = 392) channel flow, which is about three times larger than that required to sustain a turbulent 
flow.^® 

For closure at the third-moment level, turbulence transport models are important. From equation (3), 
we notice that the fourth-order moments are functions of the third-order moment through the production 
due to turbulence term [Vj^i). There is no dependence on the Reynolds stresses, except through gradients in 
the term. However, the most commonly used turbulence transport models for third-order moments are 
based on the Gaussian PDF approximation® (M) of fluctuating velocities. This model expresses the fourth- 
order moments as products of second-order moments. Figure 12 shows model comparisons with the DNS 
data. Note that in this particular case only the 122 and 222 components enter the solution, but we include 
the 111 and 112 components for completeness. The basic M model tends to over-predict the magnitude 
of the peak value for Tin and 7ii2, and the peak is predicted to be too far from the wall for 7222- This 
offset could be critical, as z;® will be overpredicted near the wall, indirectly affecting the wall shear stress, 
an important quantity in APG flows. The other models, representing corrections to the M model, do not 
improve the comparisons at all. In some cases the corrected model has no noticeable impact, while in others 
the corrected model is worse. For example, the GN model yields too small a peak magnitude for the 222 
component and the KSKl model yields too large a peak magnitude for all components. 

Fourth-order turbulence transport terms may serve as either a source or sink (hgures 7, 8). The JDJ model 
predicts the profiles reasonably well as seen in figure 13, with some errors near the local minima/maxima. 
The nearest-wall peaks in T 2222 and T 1222 are offset away from the wall. The magnitudes of the second 
(larger) peaks are overpredicted for all components except ? 2222 - The structural turbulence model (JCRl) 
presented in equation (19) for the cross-correlations uses the JDJ model to predict autocorrelations u® and 
n®. Predictions using this model are comparable to those of JDJ, with an improved prediction of the 1111 
component and under-prediction of 1222 component. Note that JCRl does not provide a prediction for 
72222- 

Third-order velocity-pressure-gradient terms are sources near the wall {yw/5 < 0.2), contributing to 
the production of v'^ and —uv^ (figures 5, 6). Further away from the wall they are a mild sink. These are 
dominant terms in the buffer region. The models NT, KSK2, DLl, and JCR2 include the net effect of a linear 
term (‘slow’ response) which is a function of the third-order moment and a non-linear term (‘fast’ response) 
which is a function of third-order moment, strain, Reynolds stresses and its gradients. Comparisons are 
plotted in figure 14. Generally, all models are poor, indicating that this component of turbulence modeling 
is in need of new ideas with improved capabilities. Recall that both NT and DLl model the pressure-strain 
term, rather than velocity-pressure-gradient. Therefore both the H and </> results from the DNS are plotted 
in figure 14 so that appropriate direct comparisons can be made. It is interesting to note that the DLl model 
predicts a profile whose general shape is more representative of H 222 than (f> 222 - As the DLl model is only a 
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Figure 12. A priori model predictions of third-order turbulence transport term (Tlji) in a fully developed 
channel at Rct = 392. The terms are non-dimensionalized by The reference statistics from DNS are 

shown as a solid line. 




Vw/5 


Figure 13. A priori model predictions of fourth-order turbulence transport term (Tljim) ^ fully developed 
channel at Ret = 392. The terms are non-dimensionalized by The reference statistics from DNS are 

shown as a solid line. 
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function of moments and budget terms of the same components, the wall limiting behavior is accurate for all 
components except Ifin, which is incorrectly predicted to be 0{y). The NT model has a ‘slow’ part identical 
to that of DLl and the model asymptote as yu, — )■ 0 is similar to DL as the ‘slow’ part is of leading order 
near the wall. This model has shortcomings similar to the DLl. As mentioned earlier, the strain-dependent 
term is identical for components 111, 112, and 122 in the NT model. This indicates that the ‘slow’ term 
plays an important role distinguishing the components. Both NT and DLl severely over-predict the peaks 
and valleys for ^m. 

Among the models for II, the KSK2 model is a blended model (equation (23)), considering a ‘slow’ part 
only in the near wall region. This part is sensitized to the wall-normal vector so as to predict the correct 
asymptotic damping as — )■ 0. The JCR2 model has a ‘fast’ term that accounts for turbulent transport 

{Tiji) in addition to the production term. This model has a ‘slow’ part similar to the DLl model and hence 
has the same near- wall shortcoming for the flm component. The JCR2 approach yields good results for all 
but the 111 component, but recall that this is strictly not a model because DNS values for Tiji are used in 
its definition. The results from JCR2 indicate that if an accurate model for Tiji were included, then it would 
yield the results shown here for 11^/. 





(b) 
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- 
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Figure 14. A priori model predictions of third-order velocity— pressure-gradient bndget term (Ilij;) in a fully 
developed channel at Re-r = 392. The terms are non-dimensionalized by u^(0)/u. The reference from DNS 
is shown as , with pressnre-strain term as — e— . 



Figure 15. A priori model predictions of fonrth-order velocity— pressure-gradient budget term (Tlijlm) in a fully 
developed channel at Re-r = 392. The terms are non-dimensionalized by u^{0)/v. The reference statistics from 
DNS are shown as a solid line. 
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The fourth-order velocity-pressure-gradient term is very similar to the second-order velocity-pressure- 
gradient term, in that it is one of the dominant terms in the balance. II 2222 produces (compare figure 7 to 
figure 9 in Ref. 7), and II 1222 produces uv^ (compare figure 8 to figure 8 in Ref. 7). As shown in figure 15, the 
JCR3 and JCR4 approaches yield generally poor results compared to DNS, even with exact representations 
of Tijim as a part of their definitions. In particular, they miss the inflectional character of the 1122 and 1222 
components, and the sign of the 2222 component is incorrect. Some of the mis-prediction is likely because 
the dissipation model is not accounted for in the approaches used. Also recall that JCR3 includes only the 
‘slow’ term. So the fact that JCR3 and JCR4 predictions are comparable indicates that the ‘fast’ term 
probably does not play a dominant role in Hijim models. The wall limiting behavior of JCR3 and JCR4 
models is accurate for all components except nmi, which is incorrectly predicted to reduce at 0{y^) as 
2 /u, — > 0 (refer to the Appendix). Again, as with the third-order velocity-pressure-gradient terms, modeling 
for the fourth-order velocity-pressure-gradient terms needs improvement. 

Third-order dissipation is small relative to other budget terms, acting as a sink in the outer layer, 
y/6 > 0.05 (figures 5, 6). The DL2 model expresses as a non-linear function of Reynolds stress and 
gradient of second-order dissipation. As noted earlier, the isotropic dissipation version DL2 cannot predict 
the inflectional profile of Siji. The anisotropic-dissipation model DL3 accounts for the change in sign only 
for components £122 and £222 (figure 16), with the change in sign occurring at a wall distance larger than 
reference DNS. The DL3 model requires accounting for second-order dissipation anisotropy either through 
an algebraic or differential model for each of its Sij components. The peak magnitudes of dissipation are 
generally overpredicted by the DL2,3 models, with the largest error in £ 222 - The KSK3 model considers Siji 
to be similar to the moment it transports in the outer layer, and near the wall similarity moments are based 
on wall- normal products as in equation (28). The near- wall behavior of is accurately predicted by KSK3, 
and it predicts all the dissipation components reasonably well in general. 



Uw / ^ Vw! ^ 


Figure 16. A priori model predictions of third-order dissipation budget term (siji) in a fully developed channel 
at Re-r = 392. The terms are non-dimensionalized by The reference statistics from DNS are shown as 

a solid line. 

Fourth-order dissipation plays a dominant ‘dissipating’ role in the transport of v'^ and uv^ (figures 7, 8). 
The GGDH model of KP2 represents £ijim as a non-linear function of third-order moments and second-order 
dissipation. Due to the dependence on the triple-correlation, the model inherently produces inflectional 
profiles (figure 17), hence incorrect. The model provides a qualitative comparison to DNS only in the region 
far away from wall. Although the near-wall asymptotic behaviour of £ijim is accurate, the magnitudes of 
dissipation are largely underpredicted. The new JGR5,6 models based on scale similarity consider Eijim to be 
scale-similar to the corresponding moment it transports. The scaling variable used for JCR5 is the second- 
order timescale, while the JCR6 model uses the fourth-order timescale variables defined in Section II. Both 
the models predict the dissipation components reasonably well, with £2222 underpredicted in magnitude. 
The JCR6 model improves the near-wall prediction of dissipation terms, especially £iin and £ 1112 - Both the 
scale-similarity models predict the correct near-wall asymptotic decay of dissipation. 
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Figure 17. A priori model predictions of fourth-order dissipation budget term {sijim) in a fully developed 
channel at Ret = 392. The terms are non-dimensionalized by w^(0)/i^. The reference statistics from DNS are 
shown as a solid line. JCR5 and JCR6 model predictions overlap in all regions except near the wall. 


B. APG-strained channel 

APG straining of the channel causes the turbulence statistics to reduce in magnitude at different rates in the 
inner and outer layers^ (figures 5-8). From equations (l)-(3), we know that the statistical moments evolve 
due to the forcing through the production due to applied strain {V^). Although the straining does not 
vary with time, the forcing is time-dependent due to non-linear dependence on moments. The objective of 
a priori testing of the strained channel is to find out if the models implicitly account for the effects of APG 
forcing through their dependence on same or lower-order moments. Statistics at a strain time ^ 22 ^ = 0.365 
(sufficiently away from the start of straining so that initial transient effects are negligible) are used to test 
the models. Similar to the unstrained channel, all u and v correlation components of budget terms are 
presented. The scales of figures 18-23 are the same as the figures in section A, in order to juxtapose model 
predictions to the unstrained case and to show reduction in moment budget magnitudes. 

Third-order turbulence transport has a reduced magnitude in the near-wall region due to applied strain 
(figure 18). The 122 and 222 components respond slower to the applied strain than the 111 and 112 
components. All the models show similar predictions relative to DNS (figure 18) compared to the unstrained 
case. The waviness in KSKl predictions are caused by gradients of third-order moments from DNS, and can 
be ignored as the gradients magnify minor statistical errors. Overall, the error in modeling predictions has 
somewhat increased due to straining. 

For the fourth-order Tijim statistics (figure 19), the JDJ model predicts the reduction in magnitude of 
all the components fairly well. The structural JCRl model, which depends on the JDJ model for autocor- 
relations, predicts the components in the outer region {y/S > 0.05) better than JDJ, but under-predicts in 
magnitude the near-wall negative peak in the 1222 component. 

Third-order velocity-pressure-gradient has a reduced magnitude in the near-wall region due to straining. 
The model of DLl and the JCR2 approach both predict this reduction to some degree; however, the errors 
in the outer layer become more noticeable and the 111 component is particularly poorly predicted near the 
wall (figure 20). Models NT and DL2 provide a better prediction of II 122 in strain than in the unstrained 
case. This could be caused by a reduced contribution from the ‘fast’ part, due to lower strain. The KSK2 
model predicts almost zero 11^7 under applied strain, which is incorrect even near the channel center. 

The fourth-order velocity-pressure-gradient magnitudes are decreased by about a factor of two by the 
effects of straining. The model predictions for the strained case (figure 21) are similar in character to the 
unstrained case, but appear to be somewhat improved relative to DNS. 

Like in the unstrained channel, third-order dissipation is a small contribution to the budget in the strained 
channel case as well (figure 5). The anisotropic model DL3 predicts a reduction in magnitude of Siji] however. 
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Figure 18. A priori model predictions of third-order turbulence transport term in a strained channel at 

A 22 ^ = 0.365. The terms are non-dimensionalized by u^{0)/v. The reference statistics from DNS are shown as 
a solid line. 
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Figure 19. A priori model predictions of fourth-order turbulence transport term (Tijim) ^ strained channel 
at A 22 t = 0.365. The terms are non-dimensionalized by u^{0)/i/. The reference statistics from DNS are shown 
as a solid line. 
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Figure 20. A priori model predictions of third-order velocity— pressure-gradient budget term in a strained 

channel at ^ 22 ^ = 0.365. The terms are non-dimensionalized by u1^(0)/u. The reference ^iji from DNS is shown 
as a solid line, with pressure-strain term, (piji as — ©— . 



Figure 21. A priori model predictions of fourth-order velocity— pressure-gradient budget term (Tlijim) in a 
strained channel at A 22 ^ = 0.365. The terms are non-dimensionalized by The reference statistics from 

DNS are shown as a solid line. 
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it is insufficient to reflect the reference DNS statistics (figure 22). Near wall prediction of £222 by DL3 has 
the same magnitude as unstrained data, indicating limited response by the model. The scale-similar KSK3 
model responds adequately to the applied strain, and predicts Siji reasonably well. 

Fourth-order dissipation budget terms reduce in magnitude significantly due to the effects of applied 
straining. The JCR5,6 models predict this behavior reasonably well, indicating that Sijim remains scale- 
similar under straining conditions. As with the unstrained channel, the KP2 model largely mis-represents 
the fourth-order dissipation terms. 
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Figure 22. A priori model predictions of third-order dissipation bndget term (ciji) in a strained channel at 
A 22 t = 0.365. The terms are non-dimensionalized by u^(0)/u. The reference statistics from DNS are shown as 
a solid line. 
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Figure 23. A priori model predictions of fourth-order dissipation budget term (eijim) in a strained channel at 
A 22 ^ = 0.365. The terms are non-dimensionalized by u^(O)/^'. The reference statistics from DNS are shown as 
a solid line. JCR5 and JCR6 model predictions overlap in all regions except near the wall. 


VII. Summary 

A DNS database for third- and fourth-order moment RANS closures has been presented for unstrained 
and strained (adverse pressure gradient) channel flow. The current DNS, whose geometry has been defined to 
allow efficient testing of RANS schemes applied to APG flows, agrees well with previously published statistics 
of second-order moment data, and includes complete budget terms for the higher-order moments. A survey 
of closure models for higher-order moments was then performed. Both existing models from the literature as 
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well as a few newly introduced ideas have been described. Most of the models for T and e relate the budgets 
to their lower-order moments through a gradient diffusion hypothesis. Velocity-pressure-gradient models 
are based on analogies to the second-order velocity-pressure-gradient model of a ‘slow’ return-to-isotropy as 
well as a ‘fast’ part. 

The extensions of these modeling ideas to higher orders have then been tested via a priori studies of 
unstrained and strained channel flow, compared to the DNS. For third-order turbulent transport, the model 
of Millionshtchikov performed fairly well, although it missed some peak values. Modifications to this model 
fared no better, and in some cases were significantly worse. For fourth-order turbulent transport, both the 
model of Jovanovic et al. as well as a new model introduced here performed about equally well overall, in 
good agreement with DNS. For the third-order velocity-pressure-gradient, all models fared poorly in general; 
based on studies using an approximated balance approach, similar difficulties are expected for fourth-order 
as well. For third-order dissipation, the model of Kawamura et al. was far superior to the other models 
investigated. For fourth-order dissipation, a new model based on a similar idea (assuming similarity of 
dissipation to the fourth-order moment) yielded the best agreement with the DNS. 

It is recognized that a priori studies are often of limited practical value, and a posteriori tests of an 
entire scheme is required to form a complete picture. Nonetheless, these current tests can perhaps be used 
as a guide to modelers to help them choose which model components might be most beneficial to forming a 
complete scheme. Also, the study demonstrates the particular need for improved model representation of the 
velocity-pressure-gradient terms. Modeling the second-order velocity-pressure-gradient is already known to 
be problematic; this difficulty clearly extends to higher-order moments as well. 

Appendix A: Near-wall asymptotics 

Very close to a solid wall the spatial rates of change of instantaneous velocity in the streamwise and span- 
wise directions are small compared with the steep variation normal to the wall. The fluctuating velocity can 
be thus conveniently expanded in a Taylor series in terms of wall-normal distance y (which, for convenience, 
replaces the used in the main text) as 

u = aiy + a 2 y'^ + 03 y^ + ... (A.l) 

v = hiy + + hy^ + .... (A. 2) 

where coefficients Oi and bi are random functions of time and the streamwise (a;) and spanwise (z) directions, 
and have zero mean values. Due to incompressibility, bi = 0. For a fully developed channel flow, it can be 
shown that U ^ y. This study is similar to the analysis by Lai and So^^ for second-moment closure budgets, 
extended here to third- and fourth-order moments and their budgets. The near-wall limiting behavior of the 
moments of interest are: 


= a\y^ + . . . 

(A.3) 

r;2 = _ 

(A.4) 

uv = a\b 2 y^ -1- . . . 

(A.5) 

v?v = a1b2y'^ + . . . 

(A.6) 

uv'^ = a\b\y^ + . . . 

(A.7) 

= b^y^ -1- . . . 

(A.8) 

rj3 = afy^ -|- . . . 

(A.9) 

u^v = a\b2y^ + . . . 

(A.IO) 

vP'v'^ = alb^y^ + . . . 

(A.11) 

rtr)3 = a\b\y'^ + . . . 

(A.12) 

= b^y^ -1- . . . 

(A.13) 

= a\y^ + . . . 

(A.14) 


The near-wall behavior of the unclosed budget terms are of importance to developing accurate turbulence 
models for predicting wall-bounded flows. With the exception of velocity-pressure-gradient (Iliji and Tlijim) 
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terms, the near-wall behavior of every other term in the moment balance equation can be analyzed, as they 
are only functions of Ui and {7^, which are known. However, the near- wall asymptotic behavior of H can 
be inferred by summing all the other budget terms and looking at the remainder. From Table A.l, we 
notice that the velocity-pressure-gradient term has the same asymptotic decay as dissipation and viscous 
diffusion, except for Hm and Hmi. Hence, similar to the understanding in second-moment near- wall 
modeling, to capture the near- wall higher moments, it may be convenient to model the terms ijliji — Siji) 
and {Uiji 

m 


Table A.l. Wall-limiting behavior of budget terms in the higher-order moment equations 
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